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CHAPTER  1 
INTRODUCTION 


1.1.  General 

The  purpose  of  this  thesis  is  to  develop  a recursive  algorithm 
for  state  estimation  for  systems  with  uncertain  models.  In  many  areas 
of  electrical  engineering,  control  and  communication  for  example,  the 
problem  of  estimating  a signal  is  a very  important  one.  When  the 
system  model  is  completely  specified  it  was  shown  [1]  that  an  optimal 
solution  can  be  obtained  under  some  assumptions  and  for  various  optimality 
criteria.  However,  in  many  practical  situations  such  a complete  descrip- 
tion is  unavailable,  therefore  the  system's  model  can  only  be  specified 
in  terms  of  the  unknown  parameters.  It  was  shown  in  [2],  [3]  among  others 
that  the  use  of  an  Incorrect  system's  model  and  noise  covariances  in 
particular  have  a dramatic  effect  on  the  system  performance  and  often  will 
give  rise  to  a problem  called  "Filter  Divergence".  When  a complete  model 
of  the  system  is  unavailable  there  are  mainly  two  methods  for  estimating 
its  states:  (a)  Adaptive  Estimation,  (b)  Mimimax  Estimation.  The 
adaptive  approach  can  be  roughly  classified  into  four  categories. 

Bayesian  [4],  [5],  Maximum  Likelihood  [6];  Correlation  [6],  [7]; 
and  Covariance  Matching  [6],  [7].  In  Bayesian  Estimation  one  tries  to 
obtain  a representation  of  the  a posteriori  density  function  viz  f(X^,0|Y^) 
where  is  the  state,  @ is  the  unknown  parameter,  and  is  the  collection 
of  observations  up  to  i.  An  a priori  knowledge-  of  the  system  and  its 
statistics  is  a necessity. 

The  Maximum  Likelihood  method  tries  to  estimate  a parameter  when 
there  is  no  a priori  knowledge  of  the  parameter  available,  i.e.  one  tries 
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to  solve  the  "Likelihood  Equation". 
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The  correlation  method  had  been  used  mainly  in  time  series 
analysis  and  econometrics.  The  basic  idea  is  to  correlate  the  output  of 
the  system  either  directly  or  after  a known  linear  operation  on  it. 

Two  different  methods  can  be  developed  by  considering 

a.  The  autocorrelation  of  the  observation  [y^j. 

b.  The  autocorrelation  of  the  innovation  process. 

Usually  the  above  schemes  have  good  large  sample  properties  but  for  small 
sample  the  adaptive  schemes  yield  state  estimates  which  have  large 
error  variance . 

A minimax  solution  to  an  estimation  problem,  minimizes  a given 
cost  function  for  the  worst  case  values  of  the  unknown  parameter.  A more 
in  depth  treatment  can  be  found  in  [1],  [10],  [11].  In  general  Minimax 
schemes  yield  too  conservative  estimates  for  large  sample  and  good 
estimates  for  small  sample. 

A third  method  which  can  be  thought  of  as  an  intermediate 
between  classical  and  adaptive  estimation  is  the  "Joint  Detection  and 
Estimation"  schemes.  Several  authors  ([12],  for  example)  have  considered 
joint  detection-estimation  schemes  in  which  it  is  assumed  that  the  signal 
to  be  estimated  is  present  with  probability  P < 1.  The  system  structure 
consist  of  a detector  which  decides  whether  a signal  is  actually  present, 
and  an  estimator  to  provide  the  estimate  when  the  detector's  decision  is 
affirmative.  The  above  scheme  assumes  a priori  knowledge  of  the 
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underlying  probabilities.  An  extension  of  the  above  approach  to  the  case 
of  M possible  modes  have  been  developed  [13].  Another  approach  suggested 
Pearson  [18]:  design  the  estimator  iL  as  a Minimax  estimator  and 
perform  the  detection  via  classical  methods.  An  extension  to  [8]  was 
suggested  by  [14],  [15]  by  assigning  a cost  to  the  entire  detection  and 
estimation  scheme  and  to  minimize  it. 

The  purpose  of  this  work  is  to  develop  a recursive  algorithm 
that  will  be  the  extension  to  [14]. 


1.2.  Problem  Statement 

Consider  the  discrete-time  linear  system  satisfying  a state 
equation  of  the  form: 

Xk+1  = 0 (k+1’k)  \ + V V k = 0,1,...  (1.2) 

and  observed  in  additive  noise 


yk+l  Sc+l^+l  + vk+l’ 


k - 0,1,... 


(1-3) 


where 

X^,  y^  are  n,  m dimensional  vectors  respectively.  The  vectors  w^,  v^, 
and  p,  and  m dimensional  uncorrelated,  white  Gaussian  noise  processes 
with  zero  mean  and  covariances  and  R^,  respectively,  where  without 
of  generality  we  can  consider  only  the  problem  where  the  unknown  parameter 
appears  only  in  Q^.  The  extension  for  the  case  in  which  the  unknown 
parameter  appears  also  in  is  straight-forward.  It  is  assumed  that 
is  of  the  following  form: 

Qv  * 0Qr 


where  Qq  a known  positive  definite  matrix,  and  9 - is  an  unknown 
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parameter  which  satisfies  at  least  one  of  the  two  bounds  under  the 
following  hypotheses: 

V 9 *b 

H1:  9 s a 

It  is  desired  to  estimate  the  state  of  the  system  based  on  the  past 
observations  [y^,  1 s i s k}  = Y^.  The  proposed  approach  attempts  to 
derive  a recursive  algorithm  for  state  estimation  in  uncertain  linear 
dynamical  systems  with  unknown  noise  parameters. 

The  proposed  scheme  attempts  to  retain  the  small  sample 
properties  of  the  standard  minimax  estimator  while  having  good  adaptive 
properties  for  large  observation  records.  For  convenience  equation  (1.2) 
may  be  rewritten  as 

*k+l  = AA  + ^ rk  wk 

where  w^  is  normalized  to  have  unit  covariance,  and  hence  r is  appro- 
priately defined. 

1 . 3.  Thesis  Outline 

In  the  following  the  proposed  approach  is  investigated.  In 
Chapter  2 we  describe  the  approach  used  to  derive  the  recursive  algorithm. 
Chapter  3 discusses  computation  procedure  and  some  relevant  computation 
aspects  of  the  problem. 
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CHAPTER  2 
PROPOSED  SCHEME 

2.1.  Proposed  Approach 

Since  our  uncertainty  on  the  unknown  parameters  appears  in  the 
form  of  bounds,  a joint  detection  and  estimation  scheme  seems  to  be  the 
most  promising.  The  joint  detection  and  estimation  approach  can  be 
described  as  follows  [14].  Let  the  measurement  space  Y^  2 R1^  be  divided 
into  two  exclusive  and  exhaustive  regions  Y^(k)  and  (k) . The  estimation 
is  to  be  performed  according  to  the  rule 

5^  = Xk(i)  if  Yk€Vk>’  i " 0.1  (2  • 1) 

A 

where  X^  is  the  estimate  of  and  X^Ci)  is  the  estimator  appropriate 
to  9 CH^,  i = 0,1  . 

A 

One  approach  suggested  by  Pearson  [8]  is  to  design  X^  to  be 
the  minimax  estimator  of  x for  9^H^  with  respect  to  some  cost  function, 
and  then  given  these  estimators  to  design  an  optimal  detector  using 
classical  detection  techniques  like  Generalized  Likelihood  Ratio. 

Another  possibility  is  to  assign  a cost  to  the  overall  joint 
detection  estimation  scheme  and  to  optimize  it.  This  approach  is  very 
appealing  since  detection  and  estimation  are  closely  related  to  each 
other.  In  fact  it  was  shown  [9],  [10]  that  both  the  causal  minimum 
variance  estimate  of  an  arbitrary  signal  process  corrupted  by  additive 
white  Caussian  noise  and  its  associated  error  covariance  matrix  can  be 
obtained  from  the  sequential  Likelihood  Ratio  by  means  of  simple  formulas. 
This  refined  approach  can  be  described  as  follows: 

Assign  an  overall  conditional  cost  denoted  by  C(9)  to  the 
overall  joint  detection  and  estimation.  One  possibility  for  a cost 


x 
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function  is  to  choose  C(9)  to  be  the  CMSE  "conditional  mean  square  error" 
defined  by 

C(9)  A E{  . i^+i  ||2 1 0}  . (2.2) 

Since  9 is  unknown  in  (2.2)  we  cannot  proceed  with  the  minimization  of 
(2.2).  Instead  a modification  to  (2.2)  which  can  be  thought  as  a 
minimax  related  approach  was  introduced  [14].  The  cost  function  will  be 
optimized  with  respect  to  the  decision  rule  (i.e.  the  estimators  appropriate 
to  9€h^,  i = 0,1;  and  to  the  detector).  We  require  that  the  estimator 

max  C(9)  , i = 0,1  (2.3) 

9 € Hi 

Since  this  problem  is  a multiple  objective  optimization  problem,  a direct 
solution  of  (2.2)  cannot  be  further  considered.  Instead  a scalar  cost 
functional  C is  to  be  minimized  and  is  defined  by  (2.4) 

C = Xn  max  C(9)  + X max  C(0)  (2.4) 

U 9 €H1  6 €h0 

which  can  be  viewed  as  a generalization  of  Magill's  performance  criterion. 
Here  X^,  X^  are  weighting  parameters  that  measure  the  relative  importance 
associated  to  errors  in  each  region  . The  analogy  can  be  drawn  as  follows 
The  case  where  the  regions  in  our  problem  formulation  are 
mutually  exclusive  and  exhaustive  is  equivalent  to  Magill's,  where  the  X^ 
represent  the  a priori  probability  that  9^  € H^.  The  objective  now  is 
to  minimize  C with  respect  to  The  cost  functional  (2.4)  was  proposed 

in  [14]  for  the  static  case.  We  consider  its  extension  to  the  dynamic 
case.  However,  both  the  decision  regions  and  the  estimators  depend  on  the 
entire  observation  record  Y^,  and  requires  increasing  computational 
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complexity  and  growing  memory  requirements  with  k.  Consequently,  a sub- 
optimal  recursive  approach  is  proposed.  Let  (2.2)  and  (2.4)  at  stage 
k+1  be  written  with  a slight  modification. 

<WS)  *E[H  Vi  - V,  <2-5) 

and 


max  C (9)  + X max  C,  (9) 
0 ^ 1 0 € Hq 


(2.6) 


where  we  assumed  without  loss  of  generality  that  \q  in  (2.4)  is  equal 
to  1. 


In  minimizing  the  expression  (2.6)  at  stage  we  assume  that 

the  earlier  decisions  remains  unchanged  i.e.  we  don't  try  to  update 

A 

the  earlier  estimates  Xj(i),  j s k,  i=0,l,  and  the  corresponding  past 

decision  regions.  We  restrict  the  problem  to  be  recursive  by  searching 

for  an  estimate  which  is  only  a function  of  the  present  observation 

^k+1  an<*  PreV:*-ous  outputs  of  the  states  of  the  states  such  as 
* 

X^(i),  i=0,l  and  the  likelihood  ratio.  As  a result  the  overall  scheme  is 
not  globally  optimal  but  stage  by  stage  optimal.  In  order  to  continue 
with  the  minimization  of  (2.6)  several  approaches  are  possible. 

1.  Proceed  with  an  exact  solution  of  (2.6) 

2.  Assume  an  estimator's  structure  and  then  proceed  with 
the  minimization. 

The  second  approach  was  chosen  mainly  because  of  the  possibility 
that  approach  (1)  can  lead  to  an  intraceable  solution  [15].  Therefore, 
assume  that  the  estimator's  structure  appropriate  to  region  9 6 H.  is 
given  by  the  following  expression,  which  is  motivated  by  the  Kalman  filter 
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Wl>  * AA(i)  + Wl)  (yk+i  - °k+iAk\(1»  1 ' O’1 

(2.7) 

We  still  haven't  said  anything  about  the  functional  relationship  between 
the  estimate  of  X denoted  X^+^_ and  the  estimator  appropriate  for 
9 € , i = 0,1  denoted  X^+^(i).  Since  we  have  shown  already  that  the 

introduction  of  a scalar  cost  functional  (2.6)  is  a generalization  of 
Magill's  performance  criterion,  we  have  decided  to  use  the  generalized 
Magi 11 's  relationship  between  the  estimate  of  X and  the  appropriate 
estimators.  Therefore,  choose  the  estimate  by 


where 


q = ^c+l(Q)  * /Vk+l(1,°-)  Xk+l(1) 

k+1  1 + A^d.O) 


Vi(1’°> 


& f(Vi|9i) 
" ^Yk+Jv 


(2.8) 


(2.9) 


is  the  likelihood-Ratio,  and  f(Y^+^|0^)  is  the  conditional  density 
function  of  the  output  given  0^.  From  the  definition  of  (2.10)  it 

is  apparent  that  the  avove  relationship  can  be  used  for  cases  where 
we  have  more  than  two  associated  regions,  therefore,  more  than  two 
estimators.  Illustrating  the  point  for  M * 3 i.e.  i * 0,1,2  we  have 


\+i 


Xk+l(°)  + Ak+l(1,0)  Xk+l(1)  +Ak+1(2,0)  Xk+1(2) 


(2.10) 


i +W1*0)  +W2’0) 

where  A ,(1,0)  and  A.  .,(2,0)  are  defined  as: 
k+l  K-ri 

1 *(*■!«!» 


Ak+l<1,0)  f(Y|S0) 


(2.11A) 
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f(Y|e2) 

Ak+1(2’0)  * f(Y le0)  (2.11B) 

In  what  follows  we  shall  denote  Ak+1(l,0)  by  Ak+1> 

2.2.  Derivation  of  the  Scheme 

Based  on  equation  (2.8)  we  can  draw  a block  diagram  of  the  system  given 
in  Fig.  1.  In  our  case  the  detector  is  a device  that  calculates  the 
relationship  given  by  (2.8)  and  also  updates  the  likelihood  ratio  A^.  The 
substitution  of  (2.8)  into  (2.6)  results  in: 


!{ll 


<W0)  + 'WW1»„2 


!•} 


+ \ max 

0 € 


E{'l  Vi  ' 


(W0) + VA  rt(1>>„2 


1 + 


\+l 


!!2I01 


(2.12) 


Utilizing  the  estimator's  structure  as  given  by  (2.7)  it  is  apparent  that 
we  can  choose  and  as  the  parameters  that  will  minimize  the  CMSE 

Ck+1  in  (2.12).  We  now  substitute  the  expressions  for  the  filtered 
estimates  Xk+^(i),  = as  defined  by  (2.7)  into  (2.14),  and  take  the 

derivatives  of  (2.12)  with  respect  to  i = 0,1  as  shown  in  detail 

in  Appendix  A to  arrive  at  the  result. 

ac 


k+1 


lE{'VA<0)*k(0)Akck+i  ' 


- <%  + \^rkVl<0)AK+i  • VH\<0>»(r(ck+l 

- Ak\(1)  wk  ~kCk+l  * <1'o+T'l)9rk”k”krkCk+l  ‘ VA(0)Vk+l 

■ eiaA(1a+i  ■ 


'oOvAV01*;'0'^  +^S«rk"A<0)AA+i  + \+A<0>AA+i 
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+ ^ ck«  \\  <°k  rk  ci;+1  + 9 ck+i  rokukuk  rk  ck« 

+^vk+i“krici;+i + ck+iAkzk(0)',k+i  ^viriv;«  + v>w 
+,,iKiiltck+iAk5ik<i»k(i>Aici+i viW0^ 

+ ck+i*A(i)»krick+i + eck+irk"k”iri;ci+i 

vk+l"krkCk+l  + Ck+lAkV1)vk+l  * ^ Ck-rtrkwkvk<.l  + vk+lvk+l 


+ x max  {e{. .. 
4 € tU 


(2.13) 


ouk+l 

~ nT  " max 
6 « «, 


W-VA<o>^‘i>Aicic+i  - VkVi)^<i)Akc; 


- \^rk“A<1)AW+i  • (l>AW+i  • 1’ivfAk\a)”kr'kck+i 

-T,2^AkV1Kr"kci+i- Vrk  Vkr'kck+i  - Vrk  v;r'kck+L 

• Vk\(0)vk+k  - Vk*k(l"w  - \^rk  Vk+l  - \^rk  Yk+1 

^Os+iW^wv™ + ^ck«rk  »A<i)ak+1  + \+A*Kcu 
+'/5ck+l\*k(0)wkr’kck+l  + 9ck+Tk  VTVk+l  + '/9vk+l”kr'kck+l 

+ Ck-lAkiEk<°>''i+l  ^k+Tk  Vktl  + Wk+l1 
+ VO Ck+lVk< 1 >*k< 1 >Vk+l+  ^ Ck+lrk”k*k(1)AkCk+l  + VlW'in 

^VA<I)';rK«  + 9 Ck+1  rk"k“k  rk  Ck+1  + vk+l  “k  rk  Ck+1 
+ Ck+lAk\(1)Vk+l  + ^ Ck+1  rk”kVk+l  + Vk+lVk+l)}} 


+ X max  [e( 


(2.14) 


i 
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Where  "Hq,  are  combinations  of  the  likelihood  ratio 


defined  as  follows: 


0 2 

«+iW 


(2.15) 


\+l 


1 <1+w2 


(2.16) 


*k+l 


(2.17) 


2 2 

<i+w 

and  TIq,  Tj  and  ^ are  their  corresponding  expected  values. 

These  equations  will  have  to  be  solved  for  i * 0,1,  to 

yield  the  desired  estimate.  In  order  to  solve  these  equations  the  expected 

values  in  (2 . 13)-(2 . 14)  have  to  be  evaluated.  Note  that  only  conditional 

expectations,  assuming  past  observations  fixed,  are  taken.  However,  the 

question  of  how  to  represent  the  likelihood-ratio  arises?  Since  we 

have  shown  before  that  the  sequential  likelihood  ratio  is  closely  related 

to  the  minimum  variance  estimate  and  its  associated  variance,  a 

sequential  likelihood  ratio  formula  given  by  [10]  was  chosen. 

Let 


Vti ' "p  (Lk+i> 


(2.18) 


then 


Lk+1  " ^ + 2 1x1  det  ^+1  “ 2 ln  det  (Ck+lPk+l|k  Ck+1  + Rk+1) 

+ 2 y;+lClyk+l  ' 2 (yk+l  " Ck+lAk\)  ,(Ck+lPk+l|kCk+l  + Rk+1) 


(yk+l  ‘ Ck+lAkXk) 


-1 


(2.19) 
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In  order  to  find  the  appropriate  expectations,  the  probability  density 
function  of  the  incremented  part  in  (2.19)  needs  to  be  derived.  For  a 
complete  analysis  of  the  derivation  of  the  density  and  the  resulting 
expected  values  involving  the  likelihood-ratio  see  Appendix  B. 

The  formula  given  in  (2.19)  has  three  parts,  the  first  depends 
only  on  the  past,  the  second  is  measurement  independent,  and  the  third 
is  measurement  dependent.  It  was  shown  (B.3)-(B.6)  that  the  likelihood 
ratio  can  be  represented  as: 

1 , f det(Ck+lPk+l|k(Q)Ck+l  + Rk+1)  I 
k+1  2 L <‘*t<ck+irk«|k<l)ci+i + W 

+ 2 (yk+l  'Ck+1AA(0))' (Ck+lPk+l|k(°)Ck+l  + Vfl^  (yk+l  ' Ck+lVV0)) 

' 2 (yk+l  ‘Ck+lAk*k(1))'(Ck+lPk+l|k(1)Ck+l  + Rk+1)  (yk+l  ‘Ck+lAk\(1)) 

(2.20) 

Equation  (2.20)  may  be  simply  expressed  as  defined  in  (B.1)-(B.7) 
which  with  (2.18)  results  in: 

\+1  - exp(vk+  | X^(0)Gk(0)Xk(0)  - | ^(DC^UX^i)}  = exp  (Vk  + Z) 

(2.21) 

Each  of  the  quadratic  forms  in  (2.20)  is  a generalized  chi-square 
with  n degrees  of  freedom,  (see  Appendix  B) . In  order  to  evaluate  the 
expected  value  of  (2.21),  or  a combination  of  it  as  appears  in  (2. 13)-(2.14) 
we  need  to  find  the  density  function  of  the  random  variable  Z which  appears 
in  the  exponent  of  (2.21)  and  has  the  form: 

Z - \ [X^(0)Gk(0)Xk(0)  -X^(l)Gk(l)Xk(l)] 


(2.22) 


i 
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The  exact  density  of  Z is  not  possible  to  express  explicitly 
hence  an  approximation  is  used  as  shown  in  Appendix  B.  The  final  asymptotic 
expression  for  the  density  of  Z is  given  by 


fi(t) 


2°1J2 

2 2 
a +o 
1 2 


2 2 
2ffl  2 
.2  2 


2° (ff1®2)nr(|) 


{•*} 


t > 0 


(2.23) 


al+a2j  2n(a1cr2)nT(|) 


til  • 


t < 0 


2 2 

where  ^ and  are  parameters  to  be  obtained  from  the  first  and  second 
moments  of  Z.  Due  to  the  lengthy  derivation  only  the  general  approach 
will  be  given  here,  and  a detailed  analysis  is  given  in  Appendix  C. 

As  shown  in  Appendix  C the  mean  and  variance  of  the  approximated  density 


(2.23)  are  given  by 


variance 


2 2 
n(<^  - crp 

4 4 

2n(c1  + a2) 


(2.24) 


(2.25) 


These  parameters  are  compared  to  the  mean  Pq  and  variance  p ^ of  Z which  can 


be  oomputed  directly  from 


= e£z]  - \ E{X^(0)Gk(0)Xk(0)  - x£(l)Gk(l)Xk(l)J 


(2.26) 


* E[Z2}  - P2  - \ E[[X^(0)Gk(0)Xk(l)  - X^(l)Gk(l)Xk(l)]2}  - P2  (2.27) 
The  exact  expressions  of  Pq  and  P^  are  given  by  (C.30)  and  (C-31). 


The  resulting  equations  become 


i 
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K 

H 

r 


2 2 
n(ai  - 


31  " 2ri(ai  + a2) 


(2.28) 

(2.29) 


where  3q,  3^  in  (2 .28) -(2 .29)  are  assumed  known.  Substituting  (2.28)  into 

2 

(2.29)  and  solving  for  a2  a ? we  arrive  at  an  equation  of  the  following  form, 
'4 

- 0 (2.30) 


+ 5 


4n 


to  be  solved  for  5* 

2 

The  solution  for  as  ln  (2.30)  may  not  be  real  which  implies 

that  our  approximation  to  the  density  of  Z is  not  adequate.  Consequently 

the  following  two  alternatives  are  considered: 

2 

(a)  The  expression  for  is  real  and  positive,  which  is  satisfied  if 

,2 


23 

3 > — 
1 n 


(2.31) 


In  this  case  the  approximation  for  the  density  is  given  by  (2.23)  with 
2 2 

parameters  a^,a2‘ 

2 

(b)  If  the  expression  for  is  complex  or  negative  the  general  model 
(2.23)  cannot  be  used  to  approximate  the  density  of  Z.  In  this  case, 
three  different  cases  are  distinguished 
(1)  If  3Q  > 0 

then  we  approximate  f (•)  by  a single  Chi -Squared  density 


a.  i 

2 1 


fz(t)  * ^72 


exp 

2“'-<r“r(S)  '2a' 


fcr}- 


t > 0 


(2.32) 


X t. 
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with  parameter 


(2)  If  PQ  < 0 

then  we  approximate  fz(*)  by  a single  Chi -Squared  density 


(2.33) 


(2.34) 


with  parameter 


(3)  If  PQ  = 0 we  shall  assume  that  P(Z*0)  * 1 Note;  that 
the  motivation  for  this  case  stems  from  the  fact  that  at  the 
initialization  stage  we  may  have 


XQ  - Xx  (2.36) 

then 

Z = j XQ(G(0)  - G(1)Xq  (2.37) 

which  in  many  practical  cases  will  either  have  G(0)  > G(l),  G(0)  < G(l), 
or  G(0)  * G(l).  It  is  unlikely  that  such  a case  will  occur  after 
the  initial  stage. 

Now  that  the  density  f (t)  has  been  determined,  we  now  proceed 

4i 

to  find  the  expected  value  of  the  functions  7]^,  7)^,  and  7^  of 
previously  defined  by  (2 . 15) -(2 . 17)  . The  expectations  of  T|^ , 7)^,  and  7^ 
are  derived  by  using  (2.21)  and  the  density  of  Z (2.23),  in  the  definition 
of  expectation. 
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>CV 


[2(cr^  -Hct^)]1172  n=0 


1 


I n(-e‘Y)  (1  + n + -±r)  , t > 0 


n-1 


2a 


n-1 


[2(a2  + a2)]n/2  n=0 


S n(-eY)  + «-l)  2 , t < 0 

=0  2af 


(2.38) 


,-Y 


zl\)  - { 


n-1 


[2 (a2  + a2)]n/2  n=0 


n(-eY)  (“^  + n - 2)  , t > 0 

2a, 


r- . 2 . 2. .n/2  . 

[2(a1+a2>]  n=0 


Yn"Vl  \2 
z n(_e  > + n;  . c < 0 
2°2 


(2.39) 


[2  (a2  + a2)]n/2  n=0  '2a 


n-1 , 


14  x i 

E n(-e  Y)  — - + n - 1 

1 


n 

2 


, t > 0 


E{T1J  - { 


[2(a‘  +a‘)]“'fc  n=0 


n 

2 


— 77?  (ri + " + v * • e < 0 (2-4o> 

2ct2 


where  the  details  may  be  found  in  Appendix  B.  Before  proceeding  to  take 
the  expected  value  of  the  terms  in  (A. 21),  (A. 23)  which  are  the  final 
expressions  for  (2.13)  and  (2.14)  a few  definitions  will  be  needed. 


£{^(0)  X^(0)}  = Pk|k(0) 
E{xk(l)  ^<1»  4 Pk|k(l) 
eCx^O)  %d>}  ^ sk 


(2.41) 

(2.42) 

(2.43) 


Using  the  above  definitions  and  the  results  obtained  for  the  expected 
value  of  the  likelihood  ratio  and  its  functions  we  proceed  to  take  the 
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the  expected  value  of  (A. 21)  and  (A. 23)  to  result  in  the  following  equations 
for  the  gains  , i = 0,1, 


° (0)  ■ ■ 1'oAkPk|k(0)AkCk+l  ‘ \AkSk  kAkCk+l  " ("o  + T,199  '"kQkrkCk+l 

,Kk+l 


+ \ICk+l[Ck+lAkPklk(0)AkCk+l  + 9 Ck+1  rk\rk°k+l  + Rk+11 

+ h^k+lVW'Vkl  + 9Ck+lWkCk«  + "k+l1 

+ xt-Vkpk|k<°>Akck+x  - VksiAicU  - (%+\)9rkQkrkck+i 

+ \Kk+l|Ck+lAkPk|k<0)AkCk+l  + 8CktlWK+l  + Rk+lI 

+ \<ll'0k+lV^+l  + 9 Ck+1  WKfl  + W 1 i <2  •44) 

° " (1)  * ‘ \AkSkAkCk+l  ” 92AkPk|kn)AkCk+l  " Pl9  Fk  VkCk+l 

5Kk+l 

- VWW+l  + ’*0Kk+lICk+lAkSkAkCk+l  + 9Ck+lWK+l 
+ w + \<«tck+iAkpkik(i>Aici,+i - 9ck«rk\ri;ck+i + 

+ ’‘t-\VkAK+l  - VkPklk<l)Ak°k+l  • YWift+l 

- VrArLcU 

+ \Kk+itck+iVkAk°k+i + 9 ViWK+i  + \+li 

+ 90Kk+l[Ck+lAkPk|  kfl)AkCk+l  + 9 Ck+1  ^Vk^+l  + \+l^ 

\> .45 ) 


After  some  algebra  these  expressions  may  be  written  compactly  in  the 
form. 


i 
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v4+i  [ai]  + ^1^+1  [biJ  * [ci5 

VkJl  ^ + Vk£  W = ^ 


(2.46) 


(2.47) 


where  the  following  notation  has  been  used. 


al  = (Gk+lAkPk|k(0)AkCk+l  + Rk+l)(1+X)  + 9lCk+irk9krkCk+l 


+ x en  c,  , ,r,  q,  r’c ' 

0 k+l  k k k k+1 


(2-48) 


bi = (ViV^+i  + W<1+X)  + 9ick+irkQkr;ck+i 


+ 90Ck+l1  kQkrkCk+l 


(2.49) 


C1  ^ “^0^  1 Pk  | k ^ ° ^ AkCk+l  ' T1lAkSkAkCk+l)(1+X)  " (T\)  + V 9 lFk V^k+l 

- x<VVeoWkck+i 


(2.50) 


a2  = (Ck+lAkSkAkCk+l  + W(1+X)  + 8lCk+1rkQkrkCk+l  + X V Ck+lWkCW 


(2.51) 


b2  (Ck+lAkPk  k(1)AkCk+l+Rk+l)  (1+X)  + Vk+^kV  kCk+l  + X90Ck+irkQkrkCk+l 


c2  ^lAkSkAkCk+l  “ ^2AkPk|k(1)AkCk+l‘^l0irkQkrkCk+l 


(2.52) 


(2.53) 


Equations  (2 .46) -(2 .47)  are  two  equations  with  two  unknown  i = 0,1 

which  may  be  solved  to  yield 


- Cfc2]  - [Cl]  fa2]}C:n0Cb2J  - \ [b^fa^1]^]}'1  (2.54) 


iv>  iVi> 


(2.55) 


( 


X 
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Now  , i - 0,1  in  (2 .46 )- (2 .47)  can  be  obtained  explicitly 

by  using  the  definitions  (2 .48)-(2 .53)  in  (2 .54)-(2 .55) . The  expressions 
found  above  allow  the  recursive  computation  of  the  filter  gains.  However, 
they  are  given  as  a function  of  the  convariances  P^|^(i)  and  S^,  so  that 
to  completely  determine  the  algorithm,  the  expressions  for  updating  the 
covariances  need  to  be  derived.  We  start  by  deriving  the  one-step 
prediction  covariances  defined  by 


VilkW-  Ef\+i|it<1>  ;k+1lk(1>’  • i"0'1 


where 


\+l|k(l)  ~ ^+1  * Ak\|k(l) 


Xk+1  " Xk+l|k('L) 


AkXk  + ^ rkWk  ‘ AkXk|k(l) 


The  substitution  of  (2.54)  into  (2.56)  yields 


(2.56) 


(2.57) 


pk+l|k(1)  • EtVk|k(1)  +-^Irk“kItVk|k<l) 

■ AkE{\ik<i>Ekik(i>lik + ww; 


AkPk|k(i)A^  + e^r^.,  i=o,i 


(2.58) 


Where  9q,  0^  are  the  values  of  9 which  maximizes  the  performance  index  in 
the  two  regions  as  specified  in  (2.6).  We  now  form  the  following  random 
variables  ^^(0)  and  Xk+^(1),  by  using  (2.7)  as  follows 
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\+i<-l)  3 \+l  " \+i(l) 

= + VoT  i kwk  - AkXk(i)  * <>^^(1)  +*/9iCk+1wk  +vfc+1] , 

i = 0,1  (2.59) 

The  substitution  of  X^^i),  i=0,l  into  the  definition  of  Sk+1  in  (2.43) 
results  after  taking  the  expectations  in  the  expressions. 

s,  ^ = a s.  a 1 + er,  q,  rv  - + K.(°h(c.  ,.a.  s.a;  + qc,  ror') 

k+1  k k k kk  k v lc+1  k+1  k+1  kkk  k+1  k k k 


+ Cwk^iti + wri + 


The  filter  covariance  error  matrix  Pk+^  is  defined  by 


(2.60) 


Pk+1  = E^*k+l*k+l^ 


(2.61) 


where 

*k+l  = \+l  ' *k+l 

= Akxk  + Je  rkwk  _ (*k+l(0)  +\+l\+l(1)) 

1 +Ak+1 

After  some  extensive  algebra  we  find  that  (2.61)  becomes 


*k+i ' tAA<0)  + rA  - <H<ck+i\V°> 

+ ^ ViVk+W  ' + ^ WA 


(2.62) 


+ Vi»/(iV 


(2.63) 


On  substitutive  (2.63)  into  (2.61),  taking  the  expectation,  simplifying, 
and  using  the  notations  "Hq,  ^ , and  T|  one  obtains 
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Pk+1  ^0Pk+l|k(0)[I"Ck+lKk+l  1 + ^k^k^k+l^+l3 
+ ^lAkSkAk  "Ck+lKk+l  )+  \Pk+l|k(1)CI"Ck+lKk+l  1 


5 .2 


+ ^o+\)i9w;^v^i)erAr^+iC  • (\+v  w^^;; 

' ^0Kk+lCk+lPk+l | k ^ ^ l^k+lCk+lAkS kAk  +%Kk+i(Ck+lPk+l  | k(0)Ck+l)Kk+l 

+ ^ l^k+lCk+lAkSkAkCk+l^k+l  ’ ^O+V^+^kV  k 

+ WlWArf-  vterf' 

+ \<Z\+A+i  + v^wSl'-  \<+icw\kk 


- v&Wk+ii  k(1)+v 

+ \Kk+lCk+lPk+l|k(I)Ck+lKk+l 

- + \<Xirk\r;cK+i^r 

+ \4ll\+AVi  + \4+i\+A?i  (2-64> 


Now,  we  have  completed  the  derivation  of  all  the  expressions 
required  for  the  computation  of  the  recursive  algorithm  and  evaluating 
its  performance.  The  detailed  computational  procedure  is  discussed  in  the 
next  chapter. 
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CHAPTER  3 

COMPUTATION  PROCEDURE 


3.1.  Introduction 

In  this  chapter,  the  computation  procedure  of  the  recursive 
algorithm  derived  in  Chapter  2,  and  other  computational  aspects  will  be 
discussed.  In  Section  3.2  the  evolution  of  the  recursive  algorithm  will 
be  considered.  Then,  the  initialization  requirements  will  be  given. 
Finally  several  simplifications  in  the  computation  are  discussed. 

3.2.  Computation  Procedure 


Step  1.  Initialization 


In  this  part  we  must  specify  initial  estimates  to:  P^|^(0), 

?k  k(l)>  sk>  ®»  ®i»  anc*  t0  t*ie  likelihood  ratio  or  more  specifically 

its  expected  value. 


Ste£_2.  Computation  of  Pfc+1|k(0)  and  pk+1 J k 

The  computation  of  the  one  step  covariance  error  matrices 
Pk+l|k(0)  and  pk+i|k(l)  is  based  on  their  equations  described  by  (2.61)- 
(2.62).  The  maximizing  values  of  9Q  and  9 , initially  will  be  the  initial 
guesses  specified  in  stage  one.  At  all  other  times  0Q  and  9^  will  be 
the  result  of  step  7 of  the  previous  iteration. 

2 2 

Step  3:  Computation  of  cr^  and 

2 2 

The  procedure  of  deriving  and  was  given  in  detail  in 
Appendix  C.  First,  the  mean  Pq  is  computed  by  using  (2.26).  It  is  a 


function  of  the  following  parameters 
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P0  = fO{Pk|k(i)’  Gk(i)’  9]  * i“°>1  (3-D 

Then  the  matrices  G^(i),  i=0,l  are  obtained  from  (2.20),  which  are  a 
function  of  the  following  parameters 

Gk(i>  = fl^Ck+l,Rk+l,Pk+l|k(l)^  1=0,1  (3,2) 


Initially  1=0>1  are  known.  In  the  computation  of  j ^ (i ) , 

i=0,l  9 is  fixed  to  its  true  value.  After  9q  and  8^  has  been  computed 
2 2 

we  can  compute  cr^,  from  (2.28)  and  (2.29). 

Step  4:  The  evaluation  of  and  the  expected  value  of  the  likelihood-ratio 

is  described  by  Equation  B.9  where  initially  we  must  use  our 
estimate  for  ^ 1*^(1)  “ 1.^(0)  w*l*-c*1  have  been  specified  in  step  1, 

with  all  the  other  parameters  known.  The  likelihood-ratio  is  computed 
via  the  recursive  relation  (B.75). 


Step  5:  Computation  of  T|q,  T)^,  T]^ 

Since  the  solution  of  and  are  given  in  terms  of  an 

infinite  series.  It  is  more  convenient  to  use  numerical  integration,  a 
32  points  Gaussian  quadrature  routine  for  example.  It  is  especially 
appealing  for  the  case  of  n * 2 since  then  the  limits  of  [0,®]  can  be 
changed  to  [0,1]  via  a change  of  variables,  otherwise  the  original 
integral 


- tv-y« 

o (l-Pe'*)2 

can  be  computed  by  placing  an  artificial  upper  limit  on  (3.3) 

•»  ux 

that  the  term  that  controls  the  convergence  of  (3.3)  is  e 


(3.3) 

by  recognizing 
Therefore 


let : 
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-u*  _ iri-a 

e =10 


(3.4) 


where  ot  is  an  arbitrarily  chosen  desired  number  sufficiently  large  in 
order  to  ensure  convergence.  Once  ot  is  chosen  x can  be  computed  to  be 
the  upper  limit  of  (3.3). 

Step  6:  Computation  of  i*0»l 

i=0,l  are  computed  using  equations  (2 .54) - (2 .55) . It 
should  be  noted  that  in  this  stage  9 is  fixed  to  its  true  value. 

Step  7:  Determination  of  max{P^+^  i’Ojl  for  0 € H 

Using  the  results  of  step  6 we  compute  the  value  of  P^+.  (i) , 
i=0,l  corresponding  to  max  9 €H^.  The  procedure  is  illustrated  for 
P^+^(0).  We  first  substitute  the  following 


k+1 


9=9. 


Pk+llk 


(1) 


into  (2.64)  and  then  find  9=8 
0=9  0 


for  which  (2.64)  is  maximized.  Fortunately,  since  P^+^(i)  is  linearly 
related  to  9 we  can  set  max  9^  to  be  at  its  boundry  values.  As  a result 
of  this  step  we  have  the  value  of  P^^i)  for  t^e  next  iteration,  and  the 
values  9q,  9^  to  be  used  in  step  2 of  the  next  iteration. 


Step  8:  The  update  of  S^+^ 

The  update  of  is  performed  via  Equation  (2.60).  Note 

that  in  this  step  9 is  fixed  to  its  true  value.  At  that  point  one 
iteration  has  been  completely  finished,  and  we  start  again  in  Step  2. 
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3.3-  On-Line  Computation  Procedure 

In  section  3.2  the  necessary  steps  for  analyzing  the  performance 
of  the  algorithm  were  described.  However,  during  the  computation  procedure 
when  measurements  are  available  the  following  changes  are  necessary. 

In  step  4 of  the  previous  section  it  is  no  longer  necessary  to 
compute  the  expected  values  of  the  likelihood-ratio  and  its  functions 
TIq,  and  . Instead  we  evaluate  using  the  sequential  likelihood 

ratio  formula  [16]  given  by  (2 .20) -(2 .21) . The  likelihood  formula  together 
with  the  derived  gains  (2.54),  (2.55)  are  then  used  in  the  expressions  of 
the  filter  (2.7)  to  obtain  the  desired  estimates.  It  should  be  noted  that 

A 

in  the  above  mode  the  observations  and  the  equations  for 

A 

and  should  be  included  appropriately. 
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CHAPTER  4 

SUMMARY  AND  CONCLUSIONS 

In  this  thesis  a joint  detection-estimation  scheme  has  been  used 
to  derive  an  on  line  recursive  algorithm  for  a class  of  systems  with 
uncertain  noise  parameters.  The  algorithm  derived  is  most  suitable  for 
a wide  class  of  systems  for  which  the  uncertainty  about  the  unknown 
parameter  appears  in  the  form  of  bounds.  It  is  also  capable  of  handling 
cases  where  the  unknown  parameters  are  time-varying  as  well  as  constant. 

The  derived  scheme  is  only  stage  by  stage  optimal  and  not  globally  optimal, 
therefore,  relieving  some  of  the  computational  requirements  such  as  memory 
requirements  and  computational  speed.  The  algorithm  is  capable  of 
approximately  computing  the  performance  of  the  system  rather  than 
simulation. 

Since  observations  are  available  in  the  on-line  filtering  mode 
we  can  compute  the  likelihood-ratio  and  its  functions  via  the  sequential 
likelihood -ratio  formula  instead  of  its  expected  value.  The  gains  of  the 
system  are  computed  directly  by  the  algorithm  and  may  in  general  depend 
on  the  observations.  It  should  be  noted  that  in  the  above  mode  the 

A 

observations  and  the  equations  for  X^^Ci)*  i=0,l  should  be  included 
appropriately.  Both  off-line  computations  of  performance  and  simulations 
have  been  made  to  demonstrate  the  information  flow  in  the  filter. 
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APPENDIX  A 


THE  DERIVATION  OF  and 


The  purpose  of  this  appendix  is  to  outline  a detailed  derivation 
of  and  up  to  the  point  where  we  need  to  take  the  expected  value. 

The  scalar  cost  functional  to  be  minimized  was  given  in  (2.12)  which 

may  be  written  as 


Vi  ’ “x„  <W9)  + x 'W9) 


(A.  1) 


0 € H 


9 € H„ 


For  convenience  and  because  of  the  symmetry  of  the  problem  only  the  first 
term  in  (A.l)  will  be  explicitly  derived.  From  (A.l)  and  (2.12)  the 
expression  for  C^+^(0)  may  be  explicitly  written  as: 

r ow°>  +WW111 

<W9>  • Ek+A+1  - 2K*1  i+k, 


\+l 

<w0)  ^k+iW1”’  (iw°>  +Ak+iVi(1))i 


<1+W 


(1+/w 


■; 


(A. 2) 


The  substitution  of  (2.7)  into  (A. 2)  results  in 


Tt+1 


<0)  - Ek+i\+r2^+i^ 

+ Ak+l/Ak*k^  +Kk+i^k+l  "Ck+lAk*k(1)))}/(1  + Ak+1) 

+ (Ak*k  (° ) + *4+1  ^ k+1  ‘ Ck+1A k*k  (° ) ) + \+l  (Ak*k ( 1 ) 

+ *4+1  ^ k+1  * Ck+lAk*k  ( l>  > )}  / ( 1+A  k+1  > ) 


(A. 3) 


1 
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We  now  need  to  take  the  derivatives  of  (A. 3)  (after  maximizing  with  respect 
to  9)  with  respect  to  and  Again  only  one  expression  will  be 

derived  due  to  the  symmetry  of  the  problem.  We  shall  use  the  following 
results  from  linear  algebra 


^ (x'Ky)  - xy' 

^ (y’Kx)  * yx’ 

^ (x'K'Ky)  = K(xy ' + yx') 


(A. 4) 
(A. 5) 

(A. 6) 


where  K is  an  nxn  matrix  and  x,y  are  n and  m column  vectors.  The  use  of 
these  relations  in  taking  the  derivatives  of  (A. 3)  with  respect  to 
and  results  in 


*W0) 

a,(0) 

9Kk+l 


Er 


f 2Xk+l  (yk+l  ' Ck-HAkXk(0)) ' , 2Ak*k(0)  (yk+l'Ck+lAk\(0) ) ' 


1 + 


V 


+1 


<1+W 


,(0) 


2Kk+l  (yk+l  "Ck+lAkXk (° ) ) (yk+l  ~Ck+lAkKk (° } } ' 

<1+Vi>2 

2 Vn*k+1 ( 1 ) (yk+l  ~Ck-HA A (° ) ) ' 

<2  +\+l>2 


} 


(A. 7) 


After  some  algebra  it  may  be  further  reduced  to 


a<W°> 


' Ett'  !Vin+Vi,+!#A(0)“C('w-cw‘kV0» 
+ 2 ViVi  < » H ‘WW1  kV0)  > ' 1 1 ' <l  + Vi>2 


(A. 8) 


The  substitution  of  (1.4)  and  (2.7)  for  i“l  into  (A. 8)  results  in: 
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9C  (0) 

30) Et'-  2(1+Ak+i)Akxk  • 2<1+Viv^  rk”k  + “A'0* 

5Kk+l 


,(0) 


+ * 2ViAkxk(1)  + 2Vi^i 


<yk+rck+i\\(1))I'  t (yk+rck+iVk(0)  > ' J I'u +'W 

In  order  to  perform  the  expectation,  the  expression  (A. 9)  is  further 
simplified  in  terms  of  the  error  estimates  defined  as 


(A. 9) 


\ “ X*.  - x. 


k "k 


(A.  10) 


and  the  system  noise  and  v . The  resulting  final  form  of  (A. 9) 


becomes 


*We>  - - - - 

ZM Wk(0)xi<0>Akck+i  - 1,iAkxk<1>lti(0)Akck+i 

A+i 

- <VV^  Vk^A+i  - ye  Akxk(0)»^c' 


k+1 


- - ci0+i'i)erkwk"krkck 


k+l 


- VkV0)vk+r  VkV1)vk+i  - <VW  rkVk 


k+l 


+ \®WkV0)Xi(0>AA+l  Ck+irk“kxk<0)AkCk+l 


+ vA<0>AA+i^9ck+iAkxk<0>''Ack+i+8ck+irk"ArK 


k+l 


+ ^k+iwAck«  + ck+iAkxk<o)vi;+i  tAwiriv»A*i’wi 


\C[ck+iAkxk<1A(0>Akc;+i+^<:k«rkVi(0)AK 


k+l 
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+ vk+A(0)Aic'k+i ck+i*kxk<i)utrtcuwck+lrk“k"icric; 


k+l 


+ Vifiti  +ck+iAA(1)vk+i  ^iVAiVw1 

(A.  11) 

where  we  have  denoted  T]  , 71  , and  T]  as  defined  in  (2 . 15) -(2 . 17) . In  a 

*W0) 

similar  fashion  the  expression  for  — 777 can  be  obtained  to  result  in 


W9) 

8<il 


• 1>lAkV°A<l)AA«  - VA(1)*k(1)AkCk+l 
- rk»A(1)AK+l  - V5  rk”k*k(1,AkCk+l 


• \ ^ W'KVU  ■ V5  Vk 


k+l 


- Wkfi'w  - Wk'ffiti 


-\\V0>vi«-YkV1)vi*i-^Vk  vi+rV5^  vi 


k+l 


^1^+1  (Gk+lAk*k (°  -*K ( 1 )AkCk+l  +y^  Ck+1  rkwk^k ( 1 )AkCk 


k+l 


+ Vk+lXk+lAkCk+l  +y^  Ck+lAk\  (° } wkT kCk+l+  0Ck+l  * kWkWkrkCk 


k+l 


+ ^ vk+i  wkr  kck+i+  ck+iA  A(0)v;+i +7§  Ck+irkVk+l+  \+lVk+l> 


^0*4+1  [Gk+lAk\(1)^k(1)AkCk+l+'^Ck+l1  kWk^k(1)AkCk 


k+l 


+ Vk+lXk(1)AkGk+l 


I 
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+ ^\+i\^<1>“krK«+9ck+irkViri;ck+i+^Vi”'rK+x 


+ ck+iAA(Uvk+i  +^ck+i  r*Vk+i+'W'fcH1 


(A. 12) 


These  are  the  expressions  to  be  used  in  (2.13)  and  (2.14) 


i 
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APPENDIX  B 

THE  EVALUATION  OF  THE  EXPECTED  VALUE  OF  THE  LIKELIHOOD-RATIO 

We  have  shown  before  that  detection  and  estimation  are  highly  related 
to  each  other  via  the  sequential  likelihood  ratio  formula.  Using  the 
formula  given  in  [16]  pp.  157-161  let  the  natural  log  of  the  likelihood 
ratio  be: 

W - <"<'W  <B-l> 


The  sequential  likelihood  ratio  is  given  by  the  following  formula 


Lk+1  * *n 


rdet(ck+ipk+ii^°>c,k+i  +IW 
ldet<ck+ipk+i|k(1>cV-  +R— > 


k+l  T k+l 


+ L, 


+ 2(yk+l  “ Ck+lAkyk(0))'(Ck+lPk+l|k(0)Ck+l  + Rk+1)  (yk+l”  Ck+lAk*k(0)) 

- I<yk+i  - ck+iAk5k<1»''ck+ipk+iik<1)ck+i  + “k+U’^k+r  ck+iAk5k<1» 


(B.2) 


Now  (B.2)  may  be  expressed  in  terms  of  the  estimation  errors 

X^(i)  = - X^(i),  i » 0,1  by  substituting  into  (B.2)  the  expressions 


(1.2)  and  (1.3)  for  X^+^  and  y^+^.  After  some  algebraic  manipulation 
we  obtain 


fdet(Ck+lPk+Hk(Q)C,k+l  + Rk+l^l 
k+l  “ W <Ck+lPk+l |k(1)C' k+l  + Rk+1 V 

+ 2(Ck+lAkyk(0)+  ^ Ck+irktuk4Vk+l)'(Ck+lPk+l 

r 

+ & ck+irk«k'mk+i) 


k<o>ck+i-*k+i> 


<ck+iAkxk<°> 
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k+irk(“k+vk+l) ' (Ck+lPk+l|  k(1)Ck+l'Hlk+l) 
k+l^k^k+l 


-1 


*(ck+lAk*k^  + ^ cu-4-iri<-U)k+vic4.i  ) (B.3) 


It  should  be  noted  that  the  first  two  terms  are  constants,  as  far  as  the 
expectation  is  concerned,  and  may  be  denoted  by: 


Yl.  = * 


pdet 

nldit 


(Ck+lPk+l|k(0)Ck+l  + Rk+l)'| 


det(Ck+lPk+l|k(1)Ck+I  + Rk+1} 


+ h 


(B.4) 


Furthermore,  the  following  notations  are  used 


Vo  ■ «Wipk+i|k(1>ciUi  + Rk+i>_1  ■ 1 ■ J’: 


(B.5) 


Xk(1>  ' <CwWl>^  Ck+irk”k«k+l>'Gk<1)<Ck+lAk\(1>+-/9  Ck+lrk"’kWk+l> 

i - 0,1  (B.5a) 

The  resulting  expression  for  the  likelihood  ratio  becomes  then 
Lk+1  ' Yk  + J X^(0)Gk(0)Xk(0)  - | X(.(l)Gk(l)Xka) 


= Yk  + Z 

where 

Z = J (X^(0)Gk(0)Xk(0)  - X-(l)Gk(l)Xk(l)) 

so  that  the  likelihood  ratio  is 

Ak+i  - exp{Yk  +Z} 


(B.6) 


(B.  7) 


Since  tuk,  vk+^  and  Xk  are  Gaussian  r.v  then  the  random  variable  Z is  a 
quadratic  expression  in  variables  which  are  conditionally  Gaussian. 
Consequently,  the  conditional  densitv  of  Zk  (over  which  the  expectation  is 
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to  be  performed)  is  a generalized  chi-square  density  function  [17], 
However,  it  is  not  possible  to  obtain  a closed  form  expression  to  the 
density,  hence  an  approximate  model  is  used.  Since  Z is  a difference  of 
two  positive  semidefinite  quadratic  forms  of  n Gaussian  variables,  we 
approximate  its  density  by  a convolution  of  two  density  functions  each 
of  which  is  chi-square  with  n degrees  of  freedom,  one  with  positive  r.v 
and  one  with  negative  r.v.  This  means  that  the  density  of  Z is  given  by 


f7(t)  = f.(z)  * f_(z) 


where 


Z'“'  ‘1 

l 1 


(B.8) 


— - i 

_ — exp(-  -~)(z)2  ‘ , z > 0 

2n/2ajr(f)  2c[ 


f1(z) 


0 , 


z < 0 


(B.9) 


2n/2aJ-(f) 


exp( — r)(-z)2 
2*2 


- - 1 

«■>  * 


f2(z) 


z < 0 


z £ 0 


(B.10) 


In  this  approximation  cr^  and  cr 2 are  parameters  to  be  determined  from  the 
moments  of  the  random  variable  Z.  The  explicit  expression  for  the 
convolution  may  be  written  as: 


fz(t) 


® - z/2 a 

5— J e 

2 (a1a2)  T (^)  max(0,t) 

2 


(z-t) 

2 n . _ 2 n , 

1 T -1  20  2 7 '1 

(z)  e (z-t)  U(z-t)dz 


t/2 at 


TJT-j  expl-z^-y- 

(2CTla2)  ^ max(0, t)  2aj(j 


2 2 n . n , 

ct.  + ct_  r -1  r -1 

:pt-z(-  7 f ))(»)  (*“t)  dz 


(B.ll) 
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two  cases  need  to  be  considered. 
For  t < 0 we  have 


t/2a: 


fz(t) 


2 , 2 

a,  + ar 


-j—  r eXpt-2(-irT^))(2) 


- -1 

,2  l 


(tof2r  r-(“)  o 


2^2 


a -i 

(z-t)t  dz, 


However,  from  [19]  the  following  integral  is  used. 


J xv  1(x+3)'Qe  ,M,Xdx  = (=j-)  2exp{^]r(v)K1  (®^> 


Vtt 


where  K is  the  modified  Bessel  function  and 

v 

| arg  p|  < rr,  Rep,  > 0,  Rev  > 0 
For  our  case  in  (B.13)  we  have 


CT1  +ct2 
- 2 2 
20^2 


Q = 1 - f 


p « -t 


so  that  the  solution  to  (B.13)  becomes 


(B.12) 

(B.13) 

(B. 14a) 

(B. 14b) 

(B. 14c) 

(B. 14d ) 
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f2(t) 


exp£t/2cr  ] 


(2axa2f  T (f)  v£ 


9 2 2 


r)] 


CT1  +a2 


a.  - a_ 
exp{_t(  2~2 — 


al  + ct2 

• r<2>*l  n 


2 2 


/ 2 2 
%1CT2 


, t < 0 


(B.15) 


Equation  (B.I5)  may  be  simplified  by  using  the  asymptotic  expansion  of 
the  Bessel  function  [20], 


K (Z)  = (— )^  eZfl  + + (y-~1)(^~9)  •} 

V ' 4z;  e U + 8Z  + 2 


21  (8Z y 


(B.H- 


2 1 

where  v fixed,  |z|  large  and  p,  = 4 v . After  substituting  v =*  ^ , and 


2 2 

0 al  + a2 

|z|  “ " 2 * 2~*2  ^ into  (B.16)  we  take  only  the  first  term  of  the 

%1ct2 

resulting  asymptotic  expansion  of  the  Bessel  function  and  substitute  it 

into  (B.15)  to  obtain  the  following: 


fz(t) 


[2(aJ  + aj)] 


- - 1 

(-t)2 * * * * *  exp^-^j] 

2^2  2 


t < 0 


(B.17a) 


Similarly  for  t > 0 we  have 


°XP(— o)  o , 

2a2  - (a‘  + at)  § - 1 

fz(t)  m~T TIT  J exp{-z—-l2  2-2  }(z)2 

2 (CT 2 ^ ^ (2)  ^ 2<^1^2 


- - 1 

2 

(z-t)  dz 


(B. 17b) 
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From  [19]  and  [2l]  the  general  solution  to  (B.12)  is  obtained  from  the 
following  integral. 

J x^  1(x-u)M'"1exp(-0x)dx  * (jp) 
u Vrr  p 

with 

[Re  p,  > 0,  Re  gu  > 0]  to  result  in 


1 

2 r (p,  )exp(~2^)K  (B.  18) 

^"2 


exp(— j) 

2 at 


fz(t) 


2n(a1a2)IVi(|) 


-ict 


,,  2 2 
(2g  la2 


. n l_ 

of  2 


al  +a2 


2 2 

n (CT1  + ct2> 

r(f)exP[-t } 

*°la2 


• K 

n 

2 


1 

2 


Ct- 


(°i 


+ a2) 


2 2 
4*la2 


} 


t > 0 


(B.19) 


Following  the  same  approach  used  for  the  case  of  t < 0,  the  substitution  of 
the  asymptotic  expression  of  the  Bessel  function  into  (B.19)  results  in 

a . ! 

V®)  * " a (t)2  exPC“%}  > c > 0 (B. 20) 

[2c\  + 4:  T *°l 

Thus  fz(t)  is  obtained  approximately  by  (B.12)  and  (B.20).  We  can  now 
proceed  to  find  the  required  expected  values  of  the  functions  of  the 
likelihood-ratio  H0>  \ and  defined  by  (2.15)-(2.17).  The  expected 
value  of  Tl0  can  be  evaluated  as  follows: 
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‘tv  * J" 


-•  (1  + eYeC)2 


0 « 

(t)dt  - J*  ( • )dt  + J ( • )dt 

-OD  0 


* \ +'Tlo 


The  first  term  T|q  may  be  written  as 

- - 1 

\ " " n / L~^—y  t 2 dt 

r 2 22  n (1  + e"e  > 2^0 

C2(orJ  ^apr  T(f)  2 


J* 

n J „ 


- - 1 

.2  1 


2 2 2 n 0 (1  + eTe  ) 

[2 (aj  +a2)]  r(f) 


(t>  -t  2 exP^  dt 


2a; 


(B.21) 


* x^V^ 


0 (1  - 3e 


dx 


(B.22) 


where  Cj  is  analogous  to  the  constant  term  in  front  of  the  integral.  In 
order  to  evaluate  (B.22)  let  |p|  < 1 and  use  a geometric  series  expansion 
to  reduce  it  to  the  following  form: 


*1  "J* 


v~l  -y.x 
x e p 

0 (1  - 0e_X)2 


dx 


j- 

n=0 


(j,(pe  X)n_1dx 


(B.23) 


which  after  interchanging  summation  and  integration  becomes 
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^n_1  J vv-l0-x(tA  + n-1)^ 
n=0  0 


\ =*  2 nPn"  J x'  ~e  ■ **  ‘Mx 


2 nen_1  (y  + n-1)^  r(v) 

n«0 


(B. 24) 


This  expression  is  used  in  (B.22)  by  recognizing  that:  0 * -e\  y = — 


2a; 


and  v * ^ to  obtain  the  final  result  for  H~ 


n A 

2 . _2.,2  n=0 


[2  (oj  + aj)] 
Similarly  for  the  second  term  of  (B.21) 


2 nC-eV1^  + n-1)  2 
*2 


(B.25) 


a - i 

.2  1 


’'o  n J*  ('t}v  t 2 exPf~ dt 

2 22  n °(1+  eve  ) 2a, 

[2(a‘  +a‘)]Z  r<f>  1 


(B.26) 


which  may  be  reduced  to 


■2y 


— / 


a . i 

.2  L 


2 o 2 n 0 (1  + e e ) 

[2<oJ  +a2>]2r(f) 


exp[-t(2  + ~)}  dt  (B.27) 


2a 


Again  recognizing  that  0 * e^,  y « 2 + -iy  and  v - £ and  applying  (B.22) 

20 1 


the  final  answer  is  obtained 
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2 n(-e^)n_1(l  + n + 2 

2 o 2 2ct  . 

[2  (a2+a2)]2  1 


(B.28) 


Equations  (B.25)  and  (B.28)  are  combined  to  yield  the  final  expression  for 
the  expectation 


[2  nC-e^W  +n  + 

2 2 2 n n*0  2ct, 

[2(a2  + a2)] r(f)  1 


2 n(-e^)n_1(n  - 1 + -±r)  2} 


(B.29) 


Following  the  same  procedure  outlined  in  (B.22)-(B.29)  we  find  that 


[2(<jJ  + a2)]2T(f) 


J-2 


te_V  J*  (t)^.  _t  2 exPC-t(-^2  + 1)}  dt 

0 (1  + e^e  V 2a? 


- - 1 


+ eY  J*  y -t  2 ^PC'tC1  + “T  >}dtl 

0 (1  + eYe  V 


2 . 2.-2  n“° 


[2(a1  + a2)]‘ 


C 2 n(-e"Y)n’l(-!?  + n-1)  2 
n-0  2a J 


+ 2 n(-eY)n‘1(-i5-  + n+l)  2 } 


(B.30) 
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2 n 

[2(c J2X  + a2)]2r(§) 


~ ~ 1 

£ J*  LtJ“v  -t  2 exp  £~2^dt 

0 (1  + e Ve  Lr  2 


J-' 


+ I ^ y -t  2 exp  Ct(“V  + 2)^dt^ 


0 (1  + eTe  ") 


[2(ct2  +a2)]2r(f) 


{ E + n-1) 


n«0 


2ct, 


• . a 

+ Z n(-e“Y)n  1(~^r  + n+1)  2 (B.31; 

n-0  2a 2 

where  the  same  procedure  of  (B.22)-(B.29)  has  been  used. 

The  expected  value  of  is  needed  mainly  for  step  4 of  the 

computation  of  the  algorithm  (see  Chapter  3,  Section  2).  We  have  denoted 


Yk  + 2 

Ak+1  " e 


(B.32) 


Yk 

where  the  term  e is  constant  as  far  as  the  expectation  is  concerned  and 
2 

e is  observation  dependent.  Therefore,  (B.32)  can  be  described  by  the 
following  recursion 


Ak+1 


AkCke 


(B. 33 ) 


where  and  are  the  first  and  second  terms  in  (2.20)  respectively. 
Using  the  results  given  by  [19]  we  obtain 
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APPENDIX  C 

THE  DERIVATION  OF  THE  MEAN  AND  VARIANCE  OF  f (t)  , 

The  purpose  of  chis  appendix  is  to  derive  the  mean  and  variance  of 

fz(t). 

The  idea  here  is  that  we  know  the  mean  and  variance  of  f^(z)  and  f2(z)  / 

which  are  chi-squared  densities.  The  density  of  the  r.v  Z (2.23)  is 
approximated  by  a convolution  of  f^(z)  and  f2(z),  which  is  a density  of  a 
r.v  Z'  given  by  a sum  of  two  independent  r.v  z ^ and  z^. 

Z'  * zL  + z2  (C.l) 

Consequently,  the  mean  of  Z'  is  the  sum  of  the  means  and  the  variance  of 

Z'  is  the  sum  of  the  variances.  We  also  can  find  the  mean  and  the 

variance  of  the  original  random  variable  Z as  defined  by  (B.6).  It  is 

important  to  note  that  the  mean  and  variance  of  Z will  be  found  as  a 

function  of  the  system  parameters.  Once  we  have  the  mean  and  variance 

of  Z'  and  Z we  compare  the  appropriate  terms  and  thus  obtain  the 
2 2 

parameters  and  cr2  of  f^(z)  and  f2(z).  Since  f^(z)  and  f2(z)  are 
chi-squared  densities,  their  mean  and  variances  are  given  by 

E(z1)  - ncr^  var^)  - 2nc^  (C.2) 

2 4 

E(z2)  - na2  var(z2)  » 2na2  (C.3) 

where  we  have  noted  that  f2(z)  is  a chi-square  on  the  negative  axis. 


I. 
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Consequently  we  have  for  Z' 

E[Z'}  - E{z1  + z2}  - na^  + (-no^) 

2 2 

- n(aj  - a‘) 

var(Z)  - var(z^)  + var(z2) 

4 4 

- 2n(ax  + a2) 


(C.4) 


We  are  now  about  to  find  the  mean  of  Z via  the  likelihood-ratio. 

It  was  shown  (B.6)-(B.7)  that  the  likelihood-ratio  can  be  written  in 
the  following  general  form. 

Ak+1  - exp{yk  + | xp0)Gk(0)xk(0)  - ■|xpl)Gk(l)xk(l>3  (C.5) 

Since  yk  is  a constant  we  need  only  find  the  expected  value  of  the 
fo1* owing  term 

E[Z}  - i E[x£(0)Gk(0)xk(0)  - x£(l)Gk(l)xk(l)}  (C.6) 


Substituting  into  xk(0)  and  xk(l)  their  corresponding  terms  as  defined  by 
(B.5a)  we  obtain 


E[Z} 


i B{(S’(0)*iC'+1  W9  ejr- 

<ck+iVk<°>  Wk 

[^k^W+i »kr'kc' 

«wwl> ck+lrk 


kc'k+W ’'i+i>«k<°» 


"k  + Vi> 


k+i  + ’W<ck<1» 


"k  + vk+i)]) 


(C.7) 
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Taking  the  expected  value  of  (C.7)  and  after  some  algebra  we  find  that 
E[Z}  is  given  by  (C.8) 

80  - e(Z]  - i tcr[Pk|k(0)AtCitlGk(0)Ck+1Ak] 

- trr«VK+iGk<0>ck+irkJ  + tr[Rk+iGk<0>] 

- «[Pk|ka)AJ[Ci+1Gk(1)Ck+1Ak]  - tr[9Qkr|,C'+1Gk(l)Ci<+1rk] 

- tr[Rk+lGk(l)]]  (C.8) 

Proceeding  to  find  the  variance  of  Z consider  the  following  term 

Bj  - var[Z}  - £ E[(x^(0)Gk(0)xk(0)  - *£(l)Ck(l)xk(l)) 

. (x^(0)Gk(0)xk(0)  - x^(l)Gk(l)xk(l)T}  - ej 

- J [E[x^(0)Gk(0)xk(0)x^(0)Gk(0)xk(0)]  - ECx^(0)Gk(0)xk(0)x^(l)Gk(l)xk(l)} 

- E{x^(l)Gk(l)xk(l)x^(0)Gk(0)xk(0)} 

+ E[x^(l)Gk(l)xk(l)x^(l)Gk(l)xk(l)}  - Bq 

- (1)  + (2)  + (3)  (C.9) 

Since  the  evaluation  of  Equation  (C.9)  is  very  lengthy,  each  term  will  be 
evaluated  individually. 

(1)  Taking  the  first  term  in  (C.9),  and  expressing  the  double  quadratic 
form  with  its  indices.  It  was  shown  [18]  that  if  x is  Gaussian  then  a 
double  quadratic  form  can  be  broken  as  follows 
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Z Gij(0)Gk^(0)[ECx.(0)xj(0))E[(xk(0)x^(0)}  -r  Ejx.CO^O^E^  (0)xy  (0)} 
i> j » 

+ E[x^(0)xk(0)}E{xi(0)x^(0)}] 

■ , * Gij(0),V0>Cpij<0>V0>  + pik<0)pjt<0)  + V°>V°>1 

1 * J » 

' , h , V0,V0)CV0,V0>  + 2pik(0)pj<<°>j  <c-10> 

performing  the  summation  only  on  the  first  term  of  (C.10) 


£ Gij(0)Gk^(0)P.j(0)Pkx(0)  = tr(G(0)P(0)tnG(0)P(0)) 

if  j 

consider  now  the  second  term  of  (C.10) 


(C.ll) 


1 t\  ,2Glj«»V'»pik“»V°> 

1 9 J 9 Xj 

summing  over  the  appropriate  indices  as  follows 


^GijPjk)Gk^PU 

C(G(0)P(0))ikGkA]Pu 


[(G(0)P(0)G(0))^]Pu 


G(0)P(0)G(0)P'(0)  = > tr[G(0)P(0)G(0)P' (0)] 


(C.  12) 


(2) 


+ E 


£ Gij(0)Gki,(l)CE(xi(0)xj(0))E(xk(l)x^(l))+E(x1(0)xk(l))E(xj(0)xA(l)) 

1»  j 

(xj(0)xk(l))E(x1(0)xi(l))] 


. ,V‘,)Vl>l:V0,Vl,  + 2Vu] 

1 9 J 9 


(C.13) 
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following  the  same  procedure  used  in  (C. 10)-(C. 12)  we  find  that  the 
solution  to  (C.13)  is  given  by 

= tr[G(0)P(0)]tr[G(l)P(l)]  + 2tr[G(0)SG(l)S ']  (C.14) 

(3) 

E Gij(l)Gk£(0){E(x.(l)Xj(l))E(xk(0)x£(0))4£(xi(l)xk(0))E(xj(l)x2(0)) 

j 

+ E(x.(l)xk(0))E(Xi(l)Xi(0))J 

« tr[G(l)P(l)]tr[G(0)P(0)]  + 2tr[G(l)S 'G(0)S]  (C.15) 

After  evaluating  terms  ( 1 ) — (3 ) we  find  that  Eq.  (C.9)  becomes 
var[2}  = tr[Gk(0)Pk(0)]tr[Gk(0)Pk(0)]  + 2tr[Gk(0)Pk<0)Gk(0)P^(0)] 

- tr[Gk(0)Pk(0)]tr[Gk(l)Pk(l)]  - 2tr[Gk(0)SkGk(l)S^] 

- tr[Gk(l)Pk(l)]tr[Gk(l)Pk(l)]  - 2tr[Gk(l)S^Gk(0)Sk] 

+ tr[Gk(l)Pk(l)]tr[Gk(l)Pk(l)]  + 2tr[Gk(l)Pk(l)Gk(l)P^(l)]  - 8q 

= (C.16) 

The  substitution  of  Pk(i),  Gk(i),  i = 0,1  and  Sk  into  (C.16)  result  in  the 
following 
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v,r{Z}  - EUifx'WA'C',^  + JX  ^r^+1  + vt+l).Gk(0) 


'(citiW0)',y’scwr‘**4'**'  > ’ hK<1>KcLi+^<rLcUi^U) 


k+11  WV^k+l'  2v  k k k+1  ’ v k^k+1  ' 'k+r 


•Gk(l)*(Ck+1Akxk(l)  +v/9  Cic+Irkujk+Vk+1)^ 

*Cj(*k<°}AkCk+l  <rkCi+l'vW),Gk(0)'(ik(0)AkCk+l  + ^Ck+lVk 
+ W - |(~xk(1)AkCk+l  U)kCk+l+vk+l),Gk(L)*(Ck+lAkxk(1) 


+ v/0  Ck+Irku;k+Vk+1^  ^ “ 0O 


(C . 17  ) 


After  some  extensive  algebra  we  arrive  to  the  final  fora  of  (C.12)  given 
in  (C.18) 

«i  - si  -e2o 

+ X l:r<Aici+iGk<0>ck+iAkpk<0,)tr<Akck+iGk<0>ck+iAkpk<0)> 

+ £f<Akck+ick<0>ck+iAkpk<0>Akck+iGk<0)ck+iAkpt<0» 

+ x "«rick«Gk<o>ck+irkGk)tr<«ri;ci;+iGk«»ck+irk<’k> 

‘V^ViV^'wwtVkWWi1*;) 

+ | tr(Gk(0)Rk+1)tr(Gk(0)Rk+1)  + \ tr<Gk<0)Rk+1Gk(0)R'+1> 

- 2 "(A'kc,k+iGk<o>ck+iAkpk<o>ct<Aici;+iGk<i>ck+iAkpk<i» 

- 2 tr<AiCi+lGk(0)Ck+lAkSkAkCk+lGk<1>Ck+lAkSi) 

- 5 “«'-kck+iGk<o>ck+irkQ,"(8r;cck+iGk<i)ck+irk\) 

- 2 tr(«2rkCi+lGkCk+irkqkricCi+lGk'I>Ck+lrk,’i> 
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' 2 “<Gk<°>‘Wtr(Gk<1>Rk+l>  - i tr(Gk(0)RkTlGk(l)R;+1) 

- I c'<AkCk-flGk<»Ck+lAkStAkCk+lGk<0»Ck+lAkSk> 

- i ct(62ric'  .G.CDC.^.r^QJ/C'  Gk(0)C.  .,rkQ;)  - i cr(Gk(l)R.  ,,Gk(l)R'  ) 


^k+i^k^^k+i1  kxk‘  k^k+rk^^k+i1  k^k^  2 
- Z er<Akck+iGk(1>ck+iAkpk(1»cr(Aick+lGk<1)ck+iAkpk(1)> 

+ 2 ”<AkCkriGk<1)Gk+lAkPk<l>AkCk+lGk(1)Ck+lAkPk<1» 

+ 4 tr(erkCk+lGk<l)Ck+lI'k<!k>,:r<erkCk+lGk(1>Ck+irk\> 

+ 5 "<«2l'^i+i=k'1)ck+irkGkrick+lGk<l)ck+irkGi> 

+ i tt(Gk(l)Rk+1)tt(Gk(l)Rk+1)  + \ «(Gk(l)Rt+1Gk(l)R’+1)  - B2 

2 2 

We  now  have  two  equations  with  two  unknowns  <j^  and  • 

n(ai  - = P0 

2n(c + <^2 ) 53  3j_ 

2 2 

The  solution  of  (C . 19 ) and  (C.20)  for  and  a2  yields 


(C . 18 1 


(C.19) 

(C.20) 


2 r/?JU 

CT’  ’ [(4^ 


2 ® 1 
a,  = C® 


^0  2 i ^0 

<2^*  + 


. (la, 


(C.21) 


(C.22) 
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